#include "TypesCL.h"


__kernel void forceCalculation(const __global real4* positions,
                               __global real4* forces,
                               __local real4* forceSums,
                               const real epsilon,
                               const real sma6,
                               const uint N)
{
  uint i = get_global_id(0);
  if (i >= N) return;

  forceSums[get_local_id(0)] =(real4)0.0;
  
  for (uint j = 0; j < N; ++j)
    {
      real4 dij = positions[j] - positions[i];
      dij.w=0.0;

      real n2sqr = dij.x*dij.x + dij.y*dij.y + dij.z*dij.z;

      real num = -1*24*epsilon*sma6*(n2sqr*n2sqr*n2sqr-2*sma6);
      real den = n2sqr*n2sqr*n2sqr*n2sqr*n2sqr*n2sqr*n2sqr;
      real expression = num/den;

      dij *= expression;

      if (i != j)
        forceSums[get_local_id(0)] += dij;
    }

  forces[get_global_id(0)] = forceSums[get_local_id(0)];
}


__kernel void updatePositions( __global real4* positions,
	      		       __global real4* old_forces,
			      const __global real4* forces,
                              const __global real4* velocities,
                              const __global real* inv_2masses,
                              const real dt,
                              const uint N)
{
    uint i = get_global_id(0);
    if (i >= N) return;

    positions[i]=positions[i]+dt*velocities[i]+forces[i]*dt*dt*inv_2masses[i];
    old_forces[i]=forces[i];
}

__kernel void updateVelocities(__global real4* velocities,
	      	               const __global real4* old_forces,
			       const __global real4* forces,
			       const __global real* inv_2masses,
                               const real dt,
                               const uint N)
{
    uint i = get_global_id(0);
    if (i >= N) return;
    velocities[i]=velocities[i]+(old_forces[i]+forces[i])*inv_2masses[i]*dt;
}